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The numerical stochastic perturbation method based on Parisi-Wu quantisation 
1 is applied to a suite of simple models to test its validity at high orders. Large 

deviations from normal distribution for the basic estimators are systematically 
found in all cases ("Pepe effect"). As a consequence one should be very careful in 
| estimating statistical errors. We present some results obtained on Weingarten's 

"pathological" model where reliable results can be obtained by an application 
of the bootstrap method. We also present some evidence that in the far less 
! trivial application to Lattice Gauge Theory a similar problem should not arise at 

rfn ! moderately high loops (up to 0(a 10 ) ). 
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1. Introduction 

In a series of papers fl], [2|, [|[ it has been shown that the technique of stochastic quanti- 
sation introduced by Parisi and Wu @] can be implemented as a practical algorithm which 
enables to reach unprecedented high orders in lattice perturbation theory (e.g. a 8 in the 
plaquette expectation value in 4-D SU(3) lattice gauge theory). Evaluating perturbative 
expansion coefficients by a Monte Carlo technique opens the back-door to a number of er- 
rors, both statistical and systematic, which should be understood and, hopefully, kept under 
control. It has been shown in Ref. || how systematic errors due to the finite volume can be 
estimated, putting a limit on the perturbative order reachable on a given lattice size. The 
aim of the present paper is to present a rather detailed analysis of statistical errors, which 
present some novel features with respect to the ordinary practice in lattice gauge Monte 
Carlo. This work was triggered by the observation made by M. Pepe of unexpected large 
discrepancies with respect to known values in the perturbative coefficients of the non-linear 
0(3)<7— model. Starting from this observation ("Pepe effect") we performed a systematic 
study of simple models with a small number of degrees of freedom trying to trace the origin 
of the discrepancies and to resolve them in a reliable way. The crucial fact that our analysis 
has uncovered is the following: the statistical nature of the processes which enter into the 
calculation of perturbative coefficients is rapidly deviating from normality as we increase 
the perturbative order, i.e. the distribution function of a typical coefficient estimator is 
strongly non-Gaussian, exhibiting a large skewness and a long tail; as a result very rare 
events give a substantial contribution to the average. A simple minded statistical analysis 
based on the assumption of normality may grossly fail to identify the confidence intervals; 
some non-parametric statistical analysis, like the bootstrap method, is necessary to assess 
the statistical error and provide reliable confidence intervals. This idea will be shown at 
work in the analysis of a lattice toy model (Weingarten's "pathological" model J7j) which 
nevertheless presents many features of interest. 

In view of this analysis one should of course worry about the results obtained in Lattice 
Gauge Theory (LGT). Our main conclusion with this respect is that one is not going to 
jeopardize the picture we drew in our previous works. As amazing as it can appear at first 
sight, the application of the method to a by far less trivial model stands on a by far firmer 
ground. We shall in fact show how the distribution function of coefficient estimators for a 
typical LGT observable does not exhibit the strongly non-Gaussian nature we find in simpler 
models. 

The content of the paper is organized as follows: in sec. 2 we recall the basis facts about 
Parisi- Wu stochastic technique applied to the numerical calculation of perturbative coeffi- 
cients in quantum field theory on the lattice (hereafter the "Parisi- WU process" or PW- 
process for short). In sec. 3 we discuss the probability distributions of the PW-process; we 
shall argue that the customary asymptotic analysis of "sum of identically distributed inde- 
pendent random variables" does not really help at high orders; in this case we shall present 
numerical evidence showing that the distribution functions still present large deviations from 
normality, in particular a whole window where the density presents a power-law rather than 
Gaussian behaviour. We present some detail on the algorithmic implementation of the PW- 
process in Sec. 4 and some numerical results in Sec. 5. A discussion of an alternate route is 
given in Sec. 6, based on Girsanov's formula. We then show (sec. 7) how a bootstrap analysis 
can be very effective in estimating confidence intervals. Finally, in sec. 8, some evidence is 
produced that both convergence time of the processes and statistical errors are under control 
for LGT. We present our conclusions in sec. 9. Some details on the bootstrap method and 
a formal analysis of convergence of perturbative correlation functions are given in appendix. 
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2. Stochastic perturbation theory 

Starting from Parisi and Wu's pioneering paper, stochastic equations have been used in 
various forms to investigate quantum field-theoretical models, both perturbatively and non- 
perturbatively. In particular a Langevin approach can be used as a proposal step subject to a 
Metropolis check to implement a non-perturbative MonteCarlo for Lattice Gauge Theories. 
We are concerned here with another approach which makes use of the Langevin equation 
to calculate weak coupling expansions. The idea [l], @ is very simple: we start from the 
Langevin equation (let us focus our attention on a simple scalar field ip) 

dip(x,t) dS 

where rj(x, t) is the standard white-noise generalized process; assuming that the action S is 
splitted into a free So and an interaction part g Si, we expand the process ip into powers in 
the coupling constant 

oo 

(2) (p(x,t) = ^2g n (p n (x,t). 

n=0 

The Langevin equation is thus translated into a hierarchical system of partial differential 
equations 

[6) at ~ d Vo (^t) v[ } 

(4) — " - = -Gq 1 yj n (aj, t) + D n (<p , . . . ,<p n -i) , for n > 1 

Go being the free propagator and D n representing source terms which can be expressed in 
terms of higher functional derivatives of the interaction term S\. Notice that the source of 
randomness is confined to the first (free-field) equation; the system can be truncated at any 
order n due to its peculiar structure (D n depends only on ip m with m < n). 

The system can be used to generate a diagrammatic expansion, as Parisi and Wu did 
for gauge theories, in terms of the free propagator G ; or, it can be studied numerically by 
simulating the white-noise process, as we have discussed in Ref.0. Any given observable 
0((p) can be expanded in powers of the coupling constant 



\n>0 / n>0 



(5) 0{ v ) EOh = > yo n fa, ...,<Pn) 

\n>0 

and its expectation value is given by 

(6) (0) = /j29 n °n)=J2 c ^9 n - 

\n>0 I n>0 

The operator O n is therefore an unbiased estimator of the n— th expansion coefficient of (O). 

It is rather straightforward to implement this idea in a practical algorithm, once the theory 
has been formulated on a space-time lattice. The application to gauge theory was presented 
in Ref. using the Langevin algorithm of Batrouni et al. Finite size errors where studied 
in a subsequent paper ||. 

The key problem we want to discuss in the present paper consists in finding a reliable way 
to estimate the statistical errors in the measure of (O n ). Admittedly, it would be desirable 
to have some analytic information on the nature of the multidimensional coupled stochastic 
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processes (^). As a substitute we choose to study some toy models where we can perform 
high statistics calculations and compare the results with the exact coefficients. 

3. Toy models. 

We have extensively studied the application of numerical stochastic perturbation theory 
to the following simple models: 

i) quartic random variable: S(<p,g) = \ip 2 + \g^ , <p 6 R. 
ii) dipole random variable: S((p,g) = [1 — cos(g ip)]/g 2 , <p G (— 7T, 7r). 
Hi) Weingarten's "pathological model" @: 

where I runs over links and p over plaquettes in a simple n-dimensional cell of a cubic 
lattice. We consider n = 2, 3, 4. 

The integration measure is Y\dip exp(-S) with the ordinary Lebesgue measure dip over all 
degrees of freedom; the expectation value of any field observable is then given by 



(O) = Z- 1 J 0[p,g] exp{-S{<p,g)} . 



The calculation of the weak coupling expansion for the typical observables can be easily 
performed, except for for the model Hi) in four dimensions, where we have been unable to 
go beyond the 5 th order. We have for instance for the "quartic" integral 

(^ 4 > = -| = l-3g-24g 2 - 297g 3 - 4896 g* - 100278 g A 

- 2450304 g 5 - 69533397 g 6 - 2247492096 g 7 - 81528066378 g 8 + 0(g 9 ) , 

for the "dipole" integral 

/ / ^ 1 9 2 9* ±7g 6 29g 8 251g w 12 , 
(cos(^x)) = l- ^r-^ + ^--T77r - T^TTTT + °(d ) , 



96 512 61440 



and finally 

dlog(Z 3 (g)) 



dg 
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for Weingarten's model. We have limited ourselves to models involving real random variables 
which make it very simple to implement the stochastic differential equation Eq.([3]). To 
reach high orders we used a symbolic language to build the right-hand side (i.e. the terms 
D n (ip , . . . , ip n _i)). In the numerical experiments n is in the range 10 to 16. 

For the sake of brevity we shall discuss our numerical experiments for the third model, 
which is the nearest to lattice gauge theory, but the qualitative aspects of the results are 
indistinguishable in the other models. 

Of course one should be aware that Weingarten's model is rather peculiar; its action 
is not bounded from below and it may make sense only in Minkowski space. From the 
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perturbative viewpoint, however, it is perfectly admissible and for imaginary g the integrals 
converge absolutely So this is an example where perturbation theory for a model living in 
Minkowski space only can be computed by the stochastic method even if the model does not 
allow for a Euclidean formulation. 

4. Algorithms 

Let us describe in some detail the numerical implementation of the stochastic differential 
equations (|3|). The free field (fo can be computed exactly (as a discrete Markov chain) since 
it is a Ornstein-Uhlenbeck process HI. For a time increment r we can write 



(7) + r) = e~ T Mt) + V 7 ! " exp(-2r) N(0, 1) 

where N(m, a) is a normal deviate with mean m and variance a. By integrating over the 
time interval [t, t + r] we get 

rt+T 

(8) cp n (t + At) = e- r (p n (t) + / exp(t' -t-r) D n (tp (t'), ... , tp n -i(t')) dt' 



which can be approximated, for instance, by the trapezoidal rule. Due to the peculiar 
nesting of the differential equations it is not necessary to perform the usual "predictor" step 
to implement the trapezoidal rule, which results in a faster algorithm. 

The bias introduced by the approximation is of order 0(t u ), with v « 2, which should be 
taken into account as it is usual with the numerical implementation of Langevin equation. 
In this case however we cannot apply a Metropolis step, as it can be done in the non- 
perturbative case. The trapezoidal rule and the absence of bias on the free field, however, 
conspire to keep the finite-r error to a small value. 

The numerical experiments have been performed by running several independent trajec- 
tories in parallel (typically 10 2 to 10 4 ); a crude estimate of autocorrelation gives a value near 
to Ti n t = \ which is the relaxation time of the Ornstein-Uhlenbeck process which drives the 
whole system. Hence we measure the observable O n every n s ki P steps with n s ki P r > 1. 

The pseudo-random numbers needed to generate the normal deviates are produced through 
two different algorithms, the lagged-Fibonacci algorithm described by Knuth || and Luscher's 
algorithm [[[(J. They do not give appreciably different results in our experiments. 



5. Results 

We report the results of experiments carried out on the model Hi) since it permits a 
graduation in dimensionality. Other models have been studied in detail and they show the 
same general pattern. The first three figures are organized as follows: the first plot is a 
log- log histogram of the raw data for a given observable (i.e. On); the third plot represents 
T Jo dt averaged over 1000 independent histories. The middle plot represents a blow up 

of the same average together with some bootstrap sample which allow to estimate confidence 
intervals, as we describe in some detail later on. 

The main observation consists in the fact that the histograms corresponding to high order 
coefficients deviate substantially from a Gaussian distribution. We do not have a complete 
characterization of the densities; however the study of these histograms sheds some light 
on the problem. We observe a rather large window where the densities are dominated by a 
power-law fall-off. Above a certain order, typically 10-14, the power is approximately x~ 2 ~ u 
with small z/, which means that there are large deviations from the average which occur as 
very rare events but contribute to the average. An example is found in Fig. [I] where the 
time average shows a sharp kink at t « 1.12e4 and a relaxation thereof. In the first picture 
the log-log-histogram is contrasted with a x~ 2 and a; -3 behaviour. Several runs on the same 
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Figure 1 . An example of a very rare event in the On history. 




model iUzd) are necessary to get a satisfactory estimate of the perturbative coefficient (here 
Cn = 388561373/4194304 « 92.64). Another experiment gives more regular histories (see 
Fig. 0). 

At order 15 we encounter a similar pattern. We observe that the large jumps at t — le4 
and t = 1.7e4 both contribute to reach a value rather close to the exact one. Were it 
not for the histograms, which, by exhibiting a large x~ 2 ~ u window, warn about very slow 
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Figure 3. A 15 experiment. 





Figure 4. A0 5 experiment. 



convergence, one would be tempted to conclude that a plateau had been reached well before 

t S3 10 4 . 

On the contrary, low order coefficients are rather close to Gaussian behaviour and the 
convergence is very fast (see Fig. ^): 

We attempted an analytical characterization of the distribution functions P n (O n < x) for 
the processes O n . The problem is rather tricky and no explicit formula was reached. One can 
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only show that if O is a generic correlation of fields, and P n (Oi, . . . , O n , t) the distribution 
function of its first n perturbative coefficients at a Langevin time t, then a limit distribution 
(for t — ► oo) always exists, and that all its moments are finite. In fact this reduces to show 
that all correlation functions 

(9) (Tl^it) 

are finite in the limit t — > oo. Details can be found in App. |A[ 

Since all moments turn out to be finite, one is tempted to use general theorems regarding 



the sum of independent identically distributed random variables |TTJ. These could in principle 
make it unnecessary to have a detailed knowledge of P n , since we are averaging over a large 
number of independent histories and the outcome should be a Gaussian distribution with 
corrections which can be parameterized and fitted to the data. According to Chebyshev and 
Petrov 

(10) f [P^\x) - G(x)}dx = ^l { 9l^l + 9lM + . . . ) 

J-J n 



where the Qi(x) are polynomials, G(u) is the normal distribution and k is the number of 
independent histories and the convergence is even uniform on x. One has for instance 

Qf\x) = A^fi^A Qi k \x) = X[ k) x(3 - x 2 ); 

Af = ^^; ^ = <(0 (fe) ) 2 >. 

Unfortunately the expansion is, at best, an asymptotic one and no useful estimate on the 
error was obtained on the basis of these formulas. This fact has triggered our interest on 
non-parametric methods of analysis which will be presented in sec. 0. 
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Figure 5. Model i) by Girsanov's formula. 
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6. GlRSANOV'S FORMULA 



The huge jumps in processes O n at high n could in principle be caused by numerical 
inaccuracy, and for some time we suspected that this could in fact happen. It was soon 
realized that the large deviations are in fact needed to reach the right average in the long 
runs, but it was nevertheless reassuring that an alternative method with zero autocorrelation 
and a totally different algorithm gave the same pattern of configurations. As a matter of 
fact the method turns out to be less efficient than the direct solution of the system given by 
Eq.(^|), but it is important in order to show that the peculiar properties of the O n histories 
are not an algorithmic artifact. 

The method applies Cameron-Martin-Girsanov formula (see e.g. [12], |13|], §3.12) 



(11) 
where 



eU(x(T)) ) = EU(y{T)) exp(£ T ) ) 



dx(t) 
dyit) 



-Ax dt + b(x, t) dt + dw(t) 
-Ay dt + dw(t) , (y(0) = x (0)) 



r 



b(y(r),r) ■ dw(r) 



\b(y(r),T)\\ 2 dr 



and E (. \ w ) denotes the average with respect to the standard Brownian process w. We apply 
the formula like in Ref . Hl2| , hence A is given by the free inverse propagator and b(x,t) is 
proportional to the coupling constant g; therefore we can explicitly expand the "exponential 
martingale" £ (t) as a power series in g and get an explicit characterization of the expansion 
coefficients: 



(12) 



(O n )= lim E (H n (h(T))I*(T)/n\ 

1 — »oo 



!) 



where 



h{T) 



h(T) 



b(y(r),T) ■ dw(r) 



T 



\b(y(r) 



T 



\dr 



and H n are Hermite polynomials. This representation of the perturbative coefficients is 
completely independent on our expansion Eq.@. It has been implemented as a numerical 
algorithm and used to estimate the perturbative expansion for model i). Each iteration 
consists in selecting a normally distributed starting point x(0) and following the evolution 
up to a time T where transient effects are sufficiently damped; since A = 1 in our model, 
T > 1 is adequate. The samples are statistically independent, modulo the quality of the 
random number generator. By monitoring the average of £(£), which should be exactly one, 
we have a handle on the accuracy of the algorithm (finite time step and statistics). 

The application of Girsanov's formula gives a cross-check on the existence of large devi- 
ations; the algorithm, however, is much more cumbersome to implement on models other 
than the simple scalar field and even in this latter case it turns out to be less efficient. 




Figure 7. Bootstrap analysis for O3. 



7. Bootstrap analysis 

The consideration of a totally trivial random variable, a high power of a simple normal 
deviate, is sufficient to convince ourselves that the occurrence of large deviation is actually 
very natural. Consider a standard Gaussian random variable x and let O = x n . The 
probability density for O contains a power-law prefactor which dominates, for high n, over 
the exponential term. It is clear, by a saddle point argument, that for n > 5 the average 
(x n ) is dominated by large deviations in the process; this appears to be a common crux 
in all models we have considered. Having established that the high order coefficients will 
suffer from large fluctuations, it is important to use reliable tools to estimate the statistical 
fluctuations. Since our estimators O n will in general have strongly non-Gaussian distribution 
functions, which moreover are only known empirically through the numerical experiments, 
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Figure 8. Bootstrap analysis for On. 



it appears a natural choice to apply some non-parametric analysis. We present here the 
results obtained by applying the bootstrap method 113] . Let us first consider experiments on 



the 3-D Weingarten's model, restricted to an elementary cube (a total of 12 link variables, 
6 plaquettes). We subdivide the data coming from several runs at the same value r = .01 
(measuring every 100 steps) in bins of N samples (say iV = 1000) and on each bin we perform 
B bootstrap replicas (in this example B = 100). Fig.[7] reports the result for O3. In this 
case the bootstrap analysis is indistinguishable from a standard gaussian analysis, which 
means that the distribution is not too far from normal. The standard gaussian analysis is 
performed by taking a weighted mean of of averages in each sample, i.e. 

E, to> tf 



(X) 



gauss 



a 



A discrepancy between these values and the bootstrap's signals a significant deviation 
from normality, like in Fig.|8|. 

The next plot refers to the estimator On. Here the deviation from normal is relevant: it is 
clearly reflected in a strong discrepancy between gaussian-like weighted mean and bootstrap 
estimate. The bootstrap apparently gives a reliable estimate for the confidence interval, as 
can be seen in the close up picture (Fig. |]). 

Finally let us give a look at the overall pattern obtained in 3-D. Fig. [TT] is remarkably 
similar to the one obtained for the plaquette expansion coefficients in Wilson's 4-D SU(3) 
Lattice Gauge Theory. 



8. Lattice Gauge Theories 

A good message from this paper is that the picture we just drew for simple models does 
not spoil the confidence we have in Numerical Stochastic Perturbation Theory for a field of 
real physical interest like Lattice Gauge Theory. As amazing as it can appear at first sight, 
this does not totally come as a surprise. In such a theory both the huge number of degrees of 
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Figure 10. Expansion coefficients for the plaquette variable in 3-D "patho- 
logical" model: *=exact, dotted line = MonteCarlo. 



freedom and their coupling make it less likely applicable the simple argument of the previous 
section referring to a high power of a normal variable. One should also keep in mind that the 
gauge degrees of freedom make it possible to keep the norms of the perturbative components 
of the field under control by providing a restoring force which does not spoil the convergence 
of the process. 

As hard as formal arguments can be and since an exact characterization of the stochastic 
processes is missing also for the simple models we considered above, we take a pragmatic 
attitude. The lesson we can learn from simple models is that monitoring frequencies his- 
tograms is a good tool in order to assess huge deviations from normality. By inspecting 
histograms from Lattice Gauge Theory simulations one can convince oneself that in this 
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Figure 11. The coefficient of g 12 for the SU(3) plaquette on a 4 4 lattice. 



case convergence is much more reliable. This of course turns out to be consistent with the 
bootstrap and standard error analysis giving the same results. For instance, the plot in 
fig. 11 refers to a £77(3) plaquette high order perturbative coefficient on a small lattice: the 
process is manifestly safe from "Pepe-effect" . 

An independent hint for the reliability of Numerical Stochastic Perturbation Theory for 
Lattice Gauge Theory comes from |TS| . A couple of new perturbative orders in the expansion 
of the plaquette in 4-D ££7(3) (which means reaching a 10 ) are shown to be fully consistent 
with both the expected renormalon behaviour and with finite size effects on top of that. A 
more organic report on the status of the art concerning the application of the stochastic 
method to LGT will be given in |TH . 



9. Conclusions 

We have performed a series of simulations to test the numerical stochastic perturbation 
algorithm previously introduced in the context of Lattice Gauge Theories. The emergence 
of large non-Gaussian fluctuations in the statistical estimators for high order coefficients 
was observed in all the simple models we considered with a similar pattern. The study of 
histograms, giving an estimate of the distribution functions for these estimators, exhibits a 
large window where exponential fall-off is masked by an approximately power-like behaviour; 
this fact entails that in any finite run there exist large deviations which are very rare but 
contribute to the average on the same foot with more frequent events. The lesson we draw 
from this experiment is twofold: it is necessary to monitor the histograms of the estimators 
in order to assess the deviation from normality; in case of large deviations, a reliable estimate 
of statistical error should be obtained by non-parametric methods, such as the bootstrap. 
These problems do not appear to plague the application of the method to Lattice Gauge 
Theories. The same viewpoint presented for simple models results in histograms hinting at 
good convergence properties, fully consistent with all our previous experience in this field. 
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Appendix A. Correlation functions at high Langevin time. 

To simplify the argument, let us consider the model (i) with a single degree of freedom 
(it may be generalized also to realistic systems). First of all we observe that any correlation 
function of free fields <po(t) may be exactly computed, in fact 

(13) (Mt) 2m ) = (2m- l)!!(l-e- 2 *) m 

converges exponentially to a limit. Then one can proceed by induction: it is sufficient to 
show that the correlation function Eq.@, which has a total perturbative order = Y2jPj> 
may be written as a finite sum of correlation functions which have a total perturbative order 
strictly lower than pt, plus correlation functions which have a total perturbative order equal 
to pt but with less free fields, plus exponentially damped terms. 

To this end it is useful to re- write the formal solution given by Eq.(|8|) in discretized time 
(t = Nt). The idea is to separate the solution into the memory terms, representing the last 
iV — 1 steps, and the most recent contribution from lower orders and noise. 

(14) (p (t) = e- T Mt-r) + V^v(t) 

ip s {t) = e-^(t-r)+r/W(t) 

where 

n— 1 i 

f U \t) = " ££ fn-l-i(t) 
i=0 j=0 

and n(t) are independent gaussian variables with mean zero and variance 2. We now should 
put Eq.([L4]) inside Eq.(^). Let's do it for a correlation function of two fields (the extension to 
a general correlation function will be considered later). For any perturbative order i,j>0 
we have: 

{^{NT) V3 {NT)) = e- 2T (<p t (NT-T)<p 3 (NT-T)} + 

+re~ T ((fi(Nr - r) / a) (iVr)) 
+re- T (^j{Nt - t) f {t) {Nr)) 
+r 2 (f^(Nr)f^(Nr)) 
Let us use the general solution of the recursive formula: 

yi = on yi-i + Pi 

which is given by 

(TV \ N / N \ 

e n & 
m=l / n=l \m=n+l / 
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In our case we simply have a m = e 2r , independent on m, while f3 m is composed of a 
linear (B m ) and a quadratic (C m ) part in r: 



re- T BmM+r'Cn 



B m {r) = (cpiiNr - r) f®(Nr)) + (^(Nr - r) f®(Nr)) 
Cm = (f^(Nr)f^(Nr)). 



As a result: 
(15) 

(16) 



(^(Nt) Vj {Nr)) = e~ 2NT ( w (0) ^,(0)) + 



+ e 



-2Nt 



N 



Y,e 2nT (re- r B n (r)+r 2 C n ) 



n=l 



We observe that B n (r) is a correlation function of total perturbative order strictly lower 
than i + j. By inductive hypothesis we know that B n (r) has a finite limit for r — >• (at 
t = Nt fixed), let's say 



(17) 



B n (t) = B t + O(t) 



and that B t has a finite limit for t — > oo. Let us parameterize the remaining dependence on 
t as: 

(18) B t = B OQ + a ie ~ 2t t p +Y, aje~ qjt 

From Eq. (P73| ) it follows that the dependence on t p is absent in correlations of free fields, but 
such correction factors can appear at higher orders. By inserting Eqs.( ]I7 , l8" ) in Eq.(|HJ) the 
geometric series can be re-summed: 



-2Nt 



(^(0)^(0)) + 



+ e 

If we take first the limit r 



-27V 7 



(rB t + 0(r 2 )) 



2t 



1-e 

at t — Nt fixed and then t 
lim (<pi(t)ipj(t)) 



fO(r) 
oo, we find 



t^oo 



1 f> 

2-000- 



Up to now, we have not considered the case of correlation functions with some free fields 
together with higher order fields. In this case the argument runs almost in the same way, 
but rf^' is substituted by y/rrj: 



(MNt)MNt)(. . .)) 



-2t 



(MNt-t)MNt-t)(...)) + 
+2v^e- r (<p (NT-T)rj(Nr) (...)) 
+r (rj(NT)rj(Nr) (...)) 



Only terms with an even number of rj at time Nt contributes. In this case the inductive 
hypothesis must be applied to correlation functions with the same total perturbative order 
but with a lower number of free fields. 
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The same argument can be applied to a general correlation function, obtaining in the limit 
t — > oo 



(19) ( n r, 



J=0 



«o(% - 1) (11^ (^o) 2 y - 

p m— 1 i / p 

Rm S S ( II V™ Vj Vm-l-i 

m=l i=0 j=0 \Z=0 

where p is the maximum perturbative order present, and <p means that a factor tp should be 
dropped in the expression. 

Induction allows us to conclude that all correlation functions reach a finite limit. 

If one goes through the previous argument by keeping the first correction in Eq.flLSD it 
is possible to show that convergence to equilibrium is dominated by an exponential factor 
which is at least e~ 2t . However, even if in correlations of free fields the dependence on t 
has exactly an exponential form (or sum of exponentials), in general (at higher orders) one 
must expect power corrections like t a e~ 2t . The precise determination of them is tedious and 
beyond our scope. 

The formula given in Eq.flT9"|) is useful also because it allows to compute not only the 
mean values of observables but also some high moments of their distributions, at least for 
the model (i). In this way we have been able to show that the large fluctuations were neither 
an artifact of the numeric simulation nor the effect of a slow convergence to equilibrium. 

Appendix B. The bootstrap algorithm 

The idea of bootstrap is very simple. Let X be a random variable and 
N— tuple of values generated by a physical process, by the stock market, by your Monte Carlo, 
whichever the case you are currently studying. In the absence of any a priori information on 
the distribution function for X, one makes the most conservative assumption, namely one 
adopts as distribution function the empirical distribution with density 

1 N 

i=i 

One then uses pb 00 t to generate other iV-tuples; all values Xj, (j = 1, . . . ,N) being equally 
probable. Having produced B such iV-tuples one can estimate various statistics of interest, 
like mean, median, quartiles etc. For instance the definition of standard deviation is simply 



(20) 



Ax 



\ 



B 



fc=i 



where (x) k is the average computed on the k-th iV-tuple and (x) is their mean value. Rec- 
ommended values for B are in the range 25 < B < 200. It is obvious that no improvement 
can be obtained on mean values; the virtue of the method should be found in the easy esti- 
mation of standard errors which do not rely on the assumption of normality. This does not 
mean of course that one should not take care of autocorrelation problems or other sources 
of statistical inaccuracies. 

The method is implemented very easily in any language. A one-liner to produce a new 
iV-tuple in matlab is the following (if the vector X contains the original data) 

XB = X(ceil(N * rand(N.l))); 
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We refer to Refs. [T7| for a thorough discussion of the method. 
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